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We show that molecular dynamics based moves in the Minima Hopping (MH) method are more efficient than 
saddle point crossing moves which select the lowest possible saddle point. For binary systems we incorporate 
identity exchange moves in a way that allows to avoid the generation of high energy configurations. Using 
this modified Minima Hopping method we reexamine the binary Lennard Jones (BLJ) benchmark system 
with up to 100 atoms and we find a large number of new putative global minima structures. 
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I. INTRODUCTION 

In a global geometry optimization one has to search 
over many local minima until one finds the global min- 
imum. Movesi^ are necessary to jump from the catch- 
ment basin of the current local minimum into another 
catchment basin. The types of moves that are chosen for 
this hopping from one catchment basin into another has 
an important effect on the efficiency of the method. The 
majority of the moves used by researches in this field, fall 
into the following five categories. 

• Random Moves: The atoms are displaced randomly 
from the current positions. If the amplitude of 
the random displacement is large enough, the sys- 
tem will relax into another local minimum during 
a standard local geometry optimization. Random 
moves are widely used both within random search 
methods^ and within basin hopping 5 . 

• Molecular dynamics moves: One of the oldest 
global optimization methods, simulated anneal- 
ing^ is based on a modified molecular dynamics 
scheme where the temperature is lowered contin- 
uously. Molecular dynamics is also used as a move 
in the Minima Hopping method 7 as well as in var- 
ious other schemes^. 

• Library based moves: For systems such as proteins, 
where one knows in advance which kind of moves 
are important, non-physical moves, based on a li- 
brary of moves can be used 10 . 

• Transition state search based moves: Starting from 
a given local minimum one searches for a neigh- 
bouring saddle point. The system is then moved 
over the saddle point and afterwards relaxed into 
the local minimum across this saddle point. This 
kind of move is used in the ART method^, which 
is primarily a method to explore potential energy 
landscapes, but which will also visit the global min- 
imum if run long enough. Saddle points are also 
determined in the TAD method^ which allows to 
determine rates of rare events. 



• Genetic algorithms^ use mutation and crossover 
operations which can also be considered as moves. 
These moves are typically strongly system depen- 
dent and the moves for clusters^ are for instance 
different to the moves used for crystal structure pre- 
diction^. 

The density of configurations is increasing exponen- 
tially with increasing energy. It is therefore important to 
search only over an energy interval above (and including) 
the global minimum which is not too large. Otherwise 
the number of configurations in this interval is so large 
that one will miss most likely the global minimum. This 
means that one should use moves that will never or only 
rarely lead into high energy structures. If moves with 
this property are not used, the majority of the configu- 
rations has to be discarded in an acceptance/rejectance 
step which can be found in most global optimization al- 
gorithms. A large rejectance ratio is however also highly 
inefficient. If one considers moves that lead from one 
minimum into a neighbouring one, it has been shown 
that molecular dynamics based moves are very efficient 16 . 
Since a molecular dynamics trajectory has a fixed and 
limited kinetic energy it follows from energy conserva- 
tion that it cannot go over barriers that are higher than 
this kinetic energy. Since the Bell Evans Polanyi princi- 
plei 7 - tells us that minima behind low energy barriers are 
on average also low in energy, molecular dynamics based 
moves do not lead into high energy structures if their ki- 
netic energy is chosen such that they can only overcome 
low energy barriers. 

In this paper we will first compare two classes of moves 
that both are able to find low energy escape paths from 
a current minimum, namely molecular dynamics based 
moves with moves that are saddle point based. In the 
second part of the paper we will discuss moves for bi- 
nary systems. We will in particular discuss under which 
circumstances moves that exchange the identity of two 
atoms are efficient. In the third part we will apply our re- 
sulting global optimization scheme to the benchmark set 
of binary Lennard Jones (BLJ) clusters^ and show that 
many global minimum structures had been overlooked in 
previous studies. The LJ potential poses the same kind 
of problems for global geometry optimization methods as 
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other more realistic potentials^ and the efficiency of an 
algorithm for an LJ system is therefore indicative of the 
success for other potentials. 



II. SADDLE POINT ESCAPE MOVES 
VERSUS MOLECULAR DYNAMICS ESCAPE 
MOVES 

We will now compare the efficiency of various moves 
within the Minima Hopping method. The original MH 
method consists of a sequence of short molecular dy- 
namics trajectories followed by local geometry optimiza- 
tions^. The molecular dynamics trajectory allows to 
cross barriers to hop from one catchment basin into an- 
other and the local geometry optimization will then bring 
the system to the bottom of the catchment basin. 

The other kind of moves are based on saddle point 
searches. Starting from a local minimum, the system 
is propagated towards a saddle point. After reaching 
such a transition state and barely crossing it, a local ge- 
ometry optimization brings the system again down to a 
new minimum. We use the dimer method^ with a few 
modifications to search the saddle points and to identify 
the transition states in this modified version of the MH 
method. 

The dimer consists of two points Ri and R2 in the 
high dimensional space in which one wants to locate the 
transition state. If the dimer midpoint is labelled as R, 
then these two points are formed according to Ri = R + 
Ai?N and R2 = R AKN, where N is the normalized 
dimer direction. The dimer method consists of basically 
two steps. In the first one, the dimer is rotated into 
a position which gives a small torsional force and which 
aligns it with an eigenmode of the Hessian. This torsional 
force is given by F 1 - = (F a -F 2 )- (Fi-F 2 |N)N, with Fi 
and F 2 being the forces at Ri and R2, respectively. Since 
the dimer midpoint remains constant during the rotation, 
the force acting on R2 can be approximated by F2 = 
2(F — Fi) in order to reduce the total number of force 
evaluations, F denoting the force at the dimer midpoint. 
In the second step the dimer is translated along F e ff, 
where this modified force is given by F e // = — (F|N)N 
if the curvature along the dimer axis is negative and by 
F eff = F - 10(F|N)N otherwise. 

Whereas the translation is always done in a straightfor- 
ward way, there exist several different ways how the rota- 
tion can be carried out. Since the rotation part requires 
considerably more force evaluations than the transla- 
tion part, choosing a good rotation method is important. 
Most methods have the tendency to align the dimer along 
the lowest curvature mode. This is due to the fact that 
only the lowest curvature mode is a minimum with re- 
spect to the curvature; all other low curvature modes 
are saddle points. Whereas this circumstance does not 
cause any problem if one is interested in finding only one 
saddle point leading out of a given minimum, it is a se- 
vere shortcoming in our case. Due to this behaviour it 



is very likely that the dimer aligns itself with the lowest 
curvature mode after a few iterations even if it was ini- 
tially aligned along another than the lowest mode and, 
as a consequence, several searches will lead to the same 
saddle point, even if they were originally started in dif- 
ferent orthogonal directions. This is a serious problem 
since we arc interested in finding many different saddle 
points leading out of a given minimum and need therefore 
a rotation method which keeps the dimer on the initially 
selected mode, thus leading to distinct saddle points. 

One possibility that was already mentioned in the orig- 
inal paper 20 is to impose the restriction that the dimer 
is orthogonal to all previous dimer directions at the min- 
imum until the dimer is aligned along the lowest mode 
itself - which will happen as we approach the saddle point 
- since then there is no more tendency for the dimer to 
switch to another mode. This procedure requires however 
the knowledge of all lower lying modes if one is inter- 
ested to follow a higher mode. Since in our new method 
we are interested in finding systematically all low lying 
saddle points around a local minimum this condition is 
automatically fulfilled. However, we also develop another 
method where it is necessary to directly follow a higher 
mode and this orthogonalization procedure is as a conse- 
quence not suited, and we therefore will look for another 
way to stay on the initially selected mode. 

Tests with several rotation methods show that using 
Direct Inversion in Iterative Subspace (DIIS)2i is most 
suitable for our purpose, since DIIS has the tendency 
to catch the nearest lying stationary point, regardless of 
whether it is a minimum, maximum or saddle point. This 
allows us to stay on the initially selected mode with high 
reliability. Approximating the error vectors by — aF^, 
where a is a constant, we move Ri according to the stan- 
dard DIIS procedure with the modification that Ri has 
to be adjusted after each step to retain the fixed dimer 
separation. However, it turns out that it is for reasons of 
stability not good to stay on the initially selected mode 
at any cost, but to follow the lowest mode at some point 
instead. This can easily be achieved by abandoning the 
DIIS rotation and using the Lanczos method thenceforth. 
This switch to the Lanczos method was done as soon as 
the second derivative of the energy with respect to the 
number of iterations became negative. 

In the present implementation of our saddle point 
searches we put the focus on reliability and not on speed, 
since we are only interested in understanding the princi- 
ple of the various types of moves. Therefore about 1000 
force evaluations are required if we want to have a suc- 
cess rate of some 99 percent. Further tuning might still 
bring down the number of force evaluations, but it seems 
unlikely that it can be reduced by one order of magni- 
tude, which would be necessary to compete with molec- 
ular dynamics based moves. So it is clear that saddle 
point based moves are only of interest in practice if the 
global minimum can be found much faster with respect 
to the number of minima that have to be visited until 
the global one is found. 
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With the exception of the moves, whose details will be 
explained below, the standard MH algorithm was used, 
i.e. new minima are accepted if they are not higher than 
Ediff in energy and the value of Ediff is adjusted by a 
feedback mechanism such that on average half of the new 
configurations are accepted. 

We use the Lennard-Jones clusters LJ55 and LJ38 as 
test systems because they behave very differently. The 
L J55 is a structure seeker for which it is very easy to find 
the global minimum. LJ^s on the other hand is a two 
funnel system for which it is surprisingly difficult to find 
its global minimum in view of its small size. 100 global 
optimization runs are done in all cases to get well defined 
average values. 



A. Crossing the lowest barrier 

The first saddle point based type of move is conceptu- 
ally simple. According to the Bell Evans Polanyi princi- 
ple, the ideal move would be a move which escapes from 
the current catchment basin by going over the lowest bar- 
rier. There are however two problems using this type of 
move. In case a minimum is visited a second time, the 
sequence of minima that are visited would repeat itself ad 
infinitum and no new minima would be visited. The sys- 
tem would not be ergodic in a certain sense. This prob- 
lem can easily be overcome by a small modification of the 
method. If a minimum is visited for the first time one es- 
capes over the lowest barrier, if it is visited a second time 
one escapes over the second lowest barrier and so on. The 
second problem is that such a type of move would be in- 
credibly expensive numerically. Doing a one sided search 
for a single saddle point requires typically already a few 
hundred force evaluations and exploring more or less all 
the saddle points around a local minimum to find out 
which is the lowest one is even much more expensive. At 
this stage we are however only interested in understand- 
ing the efficiency of a certain type of move and we will 
for the moment not care about the cost of a single move. 
We will measure the efficiency of the moves by counting 
how many local minima will be visited on average in this 
modified saddle point search based version of the Minima 
Hopping algorithm before the global minimum is found 
and we will ignore the fact that the CPU time can be very 
long due to the cost of the moves. In our implementa- 
tion of this method we perform 50 saddle point searches 
starting from the current local minimum. Out of these 
we choose the one exhibiting the lowest barrier. Since 
the saddle point searches sometimes give saddle points 
which are not connected to the initial minimum (mean- 
ing that a local geometry optimization starting at that 
saddle point would not lead back to the initial minimum), 
these barriers may be meaningless. However, we do not 
care about this fact and simply choose that saddle point 
with the lowest value. 

An overview of the performance with this method, 
which we denote as lowest barrier (LB), is given in Ta- 



ble |H We compare the performance of the LB method 
to that of the standard molecular dynamics (MD) ver- 
sion and to that of the saddle point based lowest mode 
(LM) method which will be discussed in section Hi B[ As 
usual 22 we start our molecular dynamics trajectory in a 
soft direction, i.e. in a direction with low curvature in 
order to overcome a low barrier with a small number of 
molecular dynamics steps. This direction should however 
not exactly be identical to the direction of the lowest cur- 
vature, i.e. the lowest eigenvector of the Hessian matrix, 
because we would again loose ergodicity in this way. We 
need enough randomness in the initial direction of the 
velocity vector to be able to jump into different catch- 
ment basins when we escape repeatedly from a certain 
minimum. 

As one sees, the LB method is somewhat more effi- 
cient than MD in terms of the number of distinct local 
minima that are visited before finding the global mini- 
mum for LJ^g. In terms of the total number of minima 
the MD based escapes are however more efficient. This 
comes from the fact that the MD based escapes are more 
ergodic than the LB based escapes and repeated visits of 
the same minimum are therefore less likely. Fig. [T] shows 
that a MD trajectory has a large choice for crossing very 
low barriers, i.e. there are many low lying saddle point 
around a local minimum. In the case of LJ55 the MD 
based escapes are more efficient according to both crite- 
ria. The surprising result that MD is more efficient than 
LB comes from the fact that the molecular dynamics tra- 
jectory can cross several barriers whereas in the LB based 
moves one crosses by definition only one barrier. In the 
case of the L J55 cluster crossings of several barriers are 
frequently encountered since the whole energy landscape 
is strongly 'tilted' in the direction of the global minimum. 
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minima 


minima 


minima 


minima 


SP LM 


2703.6 


523.9 


415.5 


91.9 


MD 


1030.1 


297.4 


92.3 


28.4 


SP LB 


1626.1 


268.2 


584.0 


96.3 



TABLE I. Comparison of the performance of all three Min- 
ima Hopping versions for both L J38 and LJ55 clusters, "total 
minima" are the total number of visited minima per run, "dif- 
ferent minima" indicates how many among them are different. 
"SP LM" means the saddle point version that follows the low- 
est mode, "SP LB" the one that crosses the lowest barrier. 
The data is based on 100 runs for each version. 



B. Following the lowest mode 

Even though the results for the LB case are already 
discouraging, we present a second scheme which is more 
realistic since it does not require to find all the barri- 
ers around a local minimum to make a single move. In 
this scheme we exploit the fact that there is a correla- 
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tion between the curvature of the direction into which 
we start our saddle point search and the height of the 
saddle point found. This correlation is shown in Fig. Q] 
and it is the same correlation that is also exploited when 
we start our molecular dynamics trajectories in soft direc- 
tions. In this scheme, which we denote by lowest mode 
(LM), we search for the saddle point in the softest di- 
rection if the minimum is visited for the first time, in 
the second softest direction if it is visited for the sec- 
ond time and so on. Using DIIS for the dimer rotation 
ensures that these searches will lead to distinct saddle 
points with high probability. If we find a saddle point 
we will move the system over this saddle point and per- 
form a local geometry optimization. Each mode gives 
us two degenerate directions, namely the direction of the 
positive and negative eigenvector of this mode, and we 
can therefore perform two saddle point searches for each 
mode. 

The results in Table Q] show that this approach is in 
all cases much less efficient than the molecular dynam- 
ics based moves. It is due to the fact that even if we 
start our saddle point search in a soft direction we can 
frequently obtain very high saddle points, whereas in the 
molecular dynamics based moves energy conservation will 
prevent the crossing of such high barriers. The energy of 
the molecular dynamics trajectory is usually much larger 
than needed to cross a barrier, and one would therefore 
not expect that there is a correlation between the energy 
of the trajectory crossing from one catchment basin into 
another one and the height of the saddle point that con- 
nects the minima of the two catchment basins. However, 
we found that such a correlation does indeed exist, as 
shown in Fig [2] 
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FIG. 1. A plot of the barrier heights as a function of the cur- 
vature along the direction in which the saddle point search 
was started. In addition the numbers of the modes along 
which the search was started are distinguished by colours. 
The median value for each mode (median value of the curva- 
ture as well as the median value of the barrier) is plotted with 
a large dot. One can see that it is very unlikely to find high 
barriers along the softest (lowest curvature) directions. Along 
the stiffer (higher curvature) directions one finds increasingly 
higher barriers but in addition there exist also low barriers. 
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FIG. 2. The correlation between barrier heights and kinetic 
energy during a MD Minima Hopping run. The kinetic energy 
is scaled down by a factor of three for a better visualization. 
The straight lines are convolutions of the data with a Gaus- 
sian. Fhere is a clear correlation between the barrier height 
and the kinetic energy, even if the kinetic energy is much 
higher than the barriers. 



III. MOVES FOR BINARY SYSTEMS 

In the case of binary systems one might think that 
molecular dynamics based moves are inefficient since 
atoms of a certain type will only move by some slow dif- 
fusive motion to their right place. A process which can 
bring atoms potentially faster to their right place is an 
identity exchange where the identity of two atoms that 
are possibly far apart is exchanged. In a binary Lennard- 
Jones (BLJ) cluster the identities will be denoted by A 
and B. On the other hand we know that the efficiency 
of the molecular dynamics based moves in Minima Hop- 
ping is due to the fact that high energy configurations are 
rarely visited which leads to small values of E<n / / . Ta- 
ble pi] shows these values for several BLJ systems which 
have different size-mismatch values a. The E^iff value 
is always the value which gives on average an acceptance 
ratio of 0.5 in the standard MD based version of MH. 
a = Z^-B- is the size of the larger type-B atom in a BLJ- 
cluster where the smaller type-A atoms are chosen to 
have size 1. The interaction potential is then given by 



/ \ 12 













where a and /? are the types of atoms i and j, \plo a fi 
is the equilibrium pair separation and e a p is the well 
depth of the pair potential from atoms i and j. Wc 
set e AA = £ B B = cab = e = 1 and a AB = aAA ^ BB . 
With these settings, the only free parameter besides the 
number of type-A atoms N A is u which is chosen to be 
a e {1,1.05, ...,1.3}. 

Table ITTl also shows the acceptance ratio for identity ex- 
change moves followed by a local geometry optimization, 
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E diff 


acc./rej. 


1.05 


^19^26 


0.44 


1.50 


1.15 




0.66 


0.14 


1.20 


^4-33^63 


0.37 


0.03 


1.25 


^42^58 


0.46 


0.01 



TABLE II. accept ance/rejectance ratio for identity exchange 
moves with fixed Ediff. 



if the Edi ff of the MD move for the same system is used. 
One can see that these acceptance ratios get smaller and 
smaller with increasing a. This means that in most cases 
exchanging two atom types will lead to rather high en- 
ergy configurations and is hence less efficient than MD 
moves. In nature real atoms do not only differ by size 
(e.g. covalent radius) but also by their electronic prop- 
erties. Exchange moves are therefore expected to be effi- 
cient only if the atoms are very similar in every respect. 
If the atoms are very different there is actually a strong 
driving force present in the MD moves to put the differ- 
ent types of atoms at the right positions. If the global 
minimum structure is for instance a core-shell structure 
we obtain a core-shell like structure starting from a ran- 
dom position already after some 100-1000 MD moves (see 
Fig Atomic identity exchange moves are therefore not 
only not necessary, but would even be counterproductive. 
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FIG. 3. Figure shows how fast a core shell structure of a BLJ- 
cluster is formed if one starts from a random distribution of 
atoms of type A and B. The order parameter is defined as the 
fraction of atoms which are in the core. An order parameter 
of 1 corresponds thus to the global minimum structure of a 
particular system. The evolution towards the core-shell struc- 
ture is plotted only until the order parameter reaches the value 
0.95 where the core-shell structure is already clearly visible. 
The data is obtained for the system A39.B12 with a — 1.1 

For systems with two similar types of atoms (i.e. small 
a for BLJ-systems) identity exchange moves can however 
reduce the average search time for the global minimum. 
Due to the significantly lower acceptance ratio at con- 
stant E^ff, these exchange moves can not be treated 
on an equal footing as the MD moves. In particular it 



would be too expensive to do a local geometry optimiza- 
tion after each exchange move. For this reason we have 
incorporated exchange moves in the following way in our 
Minima Hopping algorithm: After each MD move we do 
a number of exchange moves which is roughly equal to 
the number of force evaluations required in the geome- 
try optimization. If the energy of the unrelaxed config- 
uration resulting from this exchange move increases by 
less than E re \ ax with respect to the original configura- 
tion it is relaxed and taken as the result of this combined 
MD/exchange move E re i ax is the energy that is on aver- 
age gained by a local geometry optimization starting from 
a relaxed configuration where the identity of two atoms 
is exchanged. Hence the relaxed energy of the exchanged 
configuration will be on average lower then the energy 
of the original configuration. In this way the exchange 
moves can help in finding the global minimum even if 
their acceptance probability is lower than the acceptance 
probability of the MD moves by a factor which is roughly 
equal to the number of force evaluations needed by the 
MD moves. 



IV. RESULTS FOR THE BINARY LJ 
BENCHMARK SYSTEMS 

Finding the global minimum for a binary Lennard- 
Jones cluster is significantly more difficult than for a 
mono-atomic Lennard- Jones cluster. In the limit where 
the two atoms become identical (but are still slightly dif- 
ferent) each configuration becomes degenerate with a de- 
generacy of ^j^-j^r 1 where Na is the number of BLJ 
atoms of type A and Nb the number BLJ atoms of type 
B. Such configurations, that can be transformed into each 
other by identity exchanges of the atom types, are called 
homotopes^. 

If the type-A and type-B atoms have a larger size- 
mismatch er, not all ^^^U^t 1 ' homotopes are stable, but 
nevertheless a smaller number of homotopes exist. The 
existence of stable homotopes increases the number of 
local minima of binary systems compared to the mono- 
atomic case. Depending on the system size and compo- 
sition this number can thus be significantly larger. The 
difficulty of the global optimization of binary Lennard- 
Jones systems is also reflected in the data of the Cam- 
bridge Cluster database^. For mono-atomic Lennard- 
Jones systems the putative global minima up to a clus- 
ter size of 1000 atoms are listed, but for BLJ-systems 
only up to 100 atoms. Since the putative global min- 
ima structures are given for 6 different size ratios a the 
database contains 600 structures. A first computation 
of the putative global minima in this database was done 
by Doye et al^. Andrea Cassioli, Marco Locatelli and 
Fabio Schoer>2£ reexamined the problem and found nearly 
100 new putative global minima for the 600 structures in 
the database. A few new structures were also found by 
Pullan 27 . In spite of the fact that several groups have al- 
ready reexamined the database we were able to find the 
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following 17 structures which are lower in energy than 
the structures listed in the database (status June 2006): 

• a = 1.30: BLJiqq, BLJgg, BLJgg, BLJgi, BLJgg 

• a = 1.25: BLJigg, BLJgg, BLJgg, BLJg-j, BLJg§ 

• a = 1.20: BLJioq, BLJg-j, BLJge, BLJg$ 

• a = 1.15: BLJ 9i , BLJ 93 , BLJ 92 . 

The global optimization runs for systems with a 
size- mismatch ratio of a > 1.2 were done with the 
standard MH algorithm^. The putative global minima 
with a = 1.15 have been found by using the above 
mentioned identity moves which turned out to be a 
powerful additional feature to the MH algorithm for the 
global optimization of binary systems with comparable 
atomic sizes. 

We did not systematically recalculate the whole Cam- 
bridge cluster database. However, a visualization of 
the putative global minimum structures provided by the 
database revealed that several structures didn't fit into 
the "series" under the same a due to too much disorder 
of the clusters or a still incomplete separation of core and 
shell. These structures were reexamined and new ener- 
getically lower structures, that frequently had also dif- 
ferent stoichiometric compositions, were found in many 
cases. 



A. New putative global minima 

We now present the structures corresponding to the 
new putative global minima. Using the classification cri- 
teria of Doye et al£^, we will assign all cluster structures 
to their structural type families and their symmetry point 
groups. 

A new class of structures which introduces a new region 
in the structural phase diagra m 24 ; 25 for large system sizes 
N > 98 and a > 1.25 is the global minimum structure of 
BLJioo,a=i.3- The polytetrahedral structure with discli- 
nation network can be classified as 4Z14 structure with 
point group symmetry Cs- Generally, Z14-atoms are 
part of a single disclination line whereas Z15 (see FigJS]) 
and Z16 atoms act as nodes connecting 3 or 4 discli- 
nation lines respectively. Disclination lines always pass 
edges with six tetrahedra around them, see FigJS] 
The 4 atoms with coordination number Z = 14 form 4 
pairwise disconnected (single) disclination lines ending 
at 4 Z = 13 shell-atoms^. This type of structure is 
also the corresponding putative global minimum struc- 
ture Of BLJgg a= i,3, B£Jl00,cr=1.25) BLJgg tCr —\.2h & n d 

BLJq8,<7=i.25- The core of these structures consists of 42 
atoms and is completely covered by the shell atoms (pure 
core-shell). Depending on the cluster size there is only 
an absence of one or two type-B shell atoms. The puta- 
tive ground state energy of 5£Jioo,ct=i.3 is -604.796307 
in common units. 




(a) (b) (c) 

FIG. 4. 4Z14 structure of BL Jioo,<t=i.3 with Cs symmetry 
and the putative ground state energy -604.796307. a) Shows 
the whole cluster, b) Disclination network embedded into the 
core, c) Disclination network 

A summary of the structural types and energies cor- 
responding to all other minima we found is given in 
table IIIII Fig. [5] shows the A^gB^g composition 
of BL Jg7 j(7 =i.25 and how the disclinations are embed- 
ded into the whole cluster. It is a polytetrahedral 
Z15 structure with C\ symmetry and the typical Z\h- 
disclination network. The putative global minimum en- 
ergy is -578.201634. 




(a) (b) (c) 



FIG. 5. Putative global minimum structure of BZ/Jg7,o-=i.25 
represented by the A^B^e, composition, a) 3 (single) discli- 
nation lines (yellow) passing edges with 6 tetrahedra around 
them (green), connected at the Z15-node (blue), b) Disclina- 
tion network coloured, c) Shows the whole cluster in corre- 
sponding colours 

Doye et alM- describe the structural motives of other 
common binary Lennard- Jones structures up to a system 
size TV = 100. 



V. CONCLUSIONS 

Molecular dynamics based moves were found to be op- 
timal in the context of global optimization. According 
to the Bell Evans Polyani principle one should escape 
over low barriers. Even though moves that escape over 
the lowest saddle point around a local minima are con- 
sequently expected to perform better, if the number of 
visited minima is taken as the success criterion, it turns 
out that molecular dynamics based moves are more ef- 
ficient for energy landscapes which have a strong funnel 
like structure, since in this case the MD trajectory can 
cross over several saddle points. For other energy land- 
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N 


a 


e 


stoichio- 


point 


structural 








metry 


group 


type 


100 


1.3 


-604.796307 


A12-B58 


c s 


4Z14 


99 


1.3 


-597.592233 


A42-B57 


c s 


4Z14 


98 


1.3 


-590.787413 


A41B57 


Ci 


Z14 


97 


1.3 


-583.871531 




Cx 


Z14 


96 


1.3 


-576.517002 


^-41-855 


d 


Z14 


100 


1.25 


-599.264624 


A42B58 


C a 


4Z14 


99 


1.25 


-592.138846 


A42-B57 


C s 


4Z14 


98 


1.25 


-584.930661 




Cs 


4Z14 


97 


1.25 


-578.201634 


^39^58 


Ci 


Z15 


96 


1.25 


-571.389275 




Ci 


Z15 


100 


1.2 


-591.768143 


^35^65 


Ci 


Z16Z15 


97 


1.2 


-571.392434 


^33^64 


C x 


Z16Z15 


96 


1.2 


-564.674461 


^33^63 


Cx 


Z16Z15 


95 


1.2 


-557.690639 


^33^62 


Ci 


Z16Z15 


94 


1.15 


-542.476905 


^31^63 


Cx 


Z16Z15 


93 


1.15 


-535.594853 


^31^62 


Cx 


Z16Z15 


92 


1.15 


-529.190149 


^28^64 


d 


Z16Z15 



TABLE III. System size N , energies e, compositions, point 
groups and structural types of the new putative global min- 
ima. 

scapes MD based escapes are about equally efficient in 
terms of the number of distinct visited minima, but are 
more efficient in terms of the total number of visited min- 
ima, i.e. in terms of the number of local geometry opti- 
mizations. In practice moves that escape over the lowest 
saddle point are too expensive since they require a com- 
plete exploration of the potential energy surface around 
the local minimum from which one wants to escape. In 
practice the only saddle point escape moves that are af- 
fordable would be moves where one searches for a single 
saddle point in a soft direction. It turns out that these 
saddle points can sometimes be very high and the ap- 
proach is therefore much less efficient that an approach 
using molecular dynamics based moves in soft directions. 
The important difference is that by energy conservation 
the MD based trajectory can not cross over energetically 
high saddle points whereas saddle point searches in soft 
directions can give very high saddle points. In addition 
MD based escapes require significantly less force evalu- 
ations (of the order of 100) than even a single saddle 
point search. The value of the parameter Ediff in the 
Minima Hopping method is a good measure of the qual- 
ity of the moves. If moves lead on average in other low 
energy configurations Edi / / will be small and one has to 
search only over low energy structures. One has there- 
fore to search only over a number of local minima which 
is much smaller than in the case where one has to search 
in a larger energy window above the global minimum. 
According to this criterion identity exchange moves are 
in general worse than MD based moves except for very 
small values of a. If identity exchange moves are how- 
ever added as some kind of post processing step to a MD 
based move without the need of an additional geometry 
optimization for each exchange trial, the efficiency of the 
Minima Hopping method can be improved. With such 



an improved version of the Minima Hopping method we 
were able to find several new global minima structures 
for binary Lennard Jones clusters with up to 100 atoms 
and size ratios of a = 1.15. For large values of a the 
ordinary Minima Hopping method without identity ex- 
changes was used and turned out to be powerful enough 
to find new global minima structures for size rations of 
a = 1.2, 1.25 and 1.3. 

Financial support from SNF and computing time from 
CSCS are acknowledged. We gratefully acknowledge ex- 
pert discussions with Riccardo Ferrando. 
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